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ABSTRACT 


Galerkln's  finite  element  method  Is  applied  to  a  two-dimensional  heat  convection- 
diffusion  problem  arising  in  the  hydrodynamic  lubrication  of  thrust  hearings  used  in 
naval  vessels.  A  parabolized  thermal  energy  equation  for  the  lubricant,  and  thermal 
diffusion  equations  for  both  bearing  pad  and  the  collar  are  treated  together,  with 
proper  juncture  conditions  on  the  Interface  boundaries.  It  has  been  known  that  a 
numerical  instability  arises  when  the  classical  Galerkln’s  method,  which  is  equivalent 
to  a  centered  difference  approximation,  is  applied  to  a  parabolic-type  partial 
differential  equation.  Probably  the  simplest  remedy  for  this  instability  is  to  use 
a  one-sided  finite  difference  formula  for  the  first  derivative  term  in  the  finite 
difference  method.  However,  in  the  present  coupled  heat  convection-diffusion  problem 
in  which  the  governing  equation  is  parabolized  in  a  subdomain  (lubricant),  uniformly 
stable  numerical  solutions  for  a  wide  range  of  the  Peclet  number  are  obtained  in  the 
numerical  test  based  on  Galerkln's  classical  finite  element  method.  In  the  present 
numerical  computations,  numerical  convergence  errors  in  several  error  norms  are 
presented  in  the  first  model  problem.  Additional  numerical  results  for  a  more 
realistic  bearing  lubrication  problem  are  presented  for  a  second  numerical  model. 

ADMINISTRATIVE  INFORMATION 

The  work  reported  here  was  in  support  of  the  Tribology  Program,  which  is  an 
interdepartmental  effort,  under  the  David  W.  Taylor  Naval  Ship  Research  and  Development 
Center's  (OTNSRDC)  Independent  Research  Program  Task  Area  ZR-000-01-01 ,  Work  Unit 
2832-121-50. 


INTRODUCTION 

The  main  objective  of  the  present  DTNSRDC  on-going  research  in  the  Tribology 


Program  at  DTNSRDC  is  to  develop  a  reliable  tool  to  predict  the  behavior  of  thrust 


bearings  used  in  naval  vessels  over  a  wide  range  of  operating  conditions,  with  regard 
to  both  hydrodynamic  and  boundary  lubrication.  The  first  step  toward  this  goal  is  to 
improve  the  hydrodynamic  lubrication  prediction.  In  this  first-step  approach,  an 
immediate  improvement  over  previous  investigations  is  a  full  coupling  of  the  thermal 
energy  equation  in  the  lubricant  and  the  heat  diffusion  in  the  surrounding  bearing 
and  the  collar.  In  this  approach,  a  temperature  variation  is  allowed  across  the 
lubricant  film  thickness  and  proper  matching  conditions  are  imposed  on  the  interface 
boundaries  between  the  lubricant  and  adjacent  solids.  In  the  present  analysis,  the 
standard  thermal  energy  equation  is  parabolized  by  assuming  that  the  ratio  of  the 
diffusion  term  to  the  convection  term  along  the  flow  direction  is  of  small  order. 

It  has  been  known  that  for  a  parabolic-type  partial  differential  equation  an 
instability  arises  when  the  classical  Galerkin  method  is  applied,  since  this  method 
is  equivalent  to  the  centered  difference  approximation.  Probably  the  most  extensively 
studied  problem  of  this  type  is  solution  of  the  well-known  boundary  layer  equations. 
In  the  boundary  layer  equations,  most  often  a  one-sided  finite  difference  formula 
is  used,  which  is  equivalent  to  choosing  the  basis  for  the  test  function  space  to  be 
different  from  that  for  the  trial  function  space  in  the  finite  element  method.  There 
are  many  reports  describing  the  introduction  of  various  weighting  functions  in  the 
innner  product  or  the  choice  of  a  test  function  basis  different  (  that  is  asymmetric) 
from  the  trial  function  basis  [1 ,2 ,3 ,4,5 ,6,7  ].  The  choice  of  a  test  function  basis 
different  from  the  trial  function  basts  plays  a  role  in  controlling  the  degree  of 
"upwinding"  to  maximize  accuracy.  Recently,  extensive  studies  have  been  made  for 
general  convection-diffusion  equations.  However,  rigorous  investigations  into  the 
control  of  the  degree  of  upwinding  are  limited  to  a  simple  one  dimensional  model 
problem  with  constant  coefficients.  For  more  general  situations,  maximization  of 
accuracy  by  controlling  the  degree  of  upwinding  is  not  straightforward.  Therefore, 
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it  is  desirable  to  have  a  simple  numerical  scheme  which  is  uniformly  stable  over 
a  wide  range  of  values  of  the  Peclet  number. 

In  this  report  numerical  results  are  obtained  by  the  classical  Galerkin  method. 

The  present  numerical  scheme  gives  uniformly  stable  results  over  a  wide  range  of 
Peclet  number  for  the  parabolized  thermal  energy  equation  in  the  lubricant.  The  present 
numerical  scheme  is  equivalent  to  the  centered  difference  approximations  for  both  the 
first  and  second  derivative  terms  in  the  original  differential  equation. 


MATHEMATICAL  MODELLING  OF  THE  PROBLEM 

The  computation  domain  consists  of  three  rectangul?  .nbdomains,  having  heights 

HI,  H2 ,  H3,  and  length  L  as  shown  in  Figure  1.  However,  '  choice  of  actual  boundary 

geometry  is  general  (  not  necessarily  rectangular  domai'  in  the  present  method.  The 
upper,  middle,  and  lower  subdomains  are  the  bearing  pad  .  lubricant,  and  the  bearing 
collar,  denoted  by  (2V,  ,  and  fl3  ,  respectively.  A  rectangular  co-ordinate  system 

is  used,  the  y-axis  pointing  upward  and  the  x-axis  pointing  toward  the  right-hand  side. 
The  origin  is  taken  at  the  mid-point  of  the  left-hand  side  vertical  boundary  of  the 
lubricant.  The  boundary  of  each  subdomain  is  as  shown  in  Figure  1, 

90;  -  T[  U  Ji  (i  =  l  and  3) 

3  Qz  —  O  rip  J  i  KJ  J  x 


(1) 


Here  the  flow  velocity  U  is  also  shown  in  Figure  1. 

In  the  present  numerical  test,  we  ignore  the  convection  term  in  the  solids,  i.e., 
ft  I  and  Q3  ,  by  assuming  they  are  stationary.  The  thermal  conductivity,  kj  (1=1,3) 
is  assumed  constant  in  each  subdomain.  In  the  subdomain  ft2  (the  lubricant),  the 
velocity  distribution  is  specifed  a  priori,  hence  the  heat  generating  source  term 
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a. 


n 


Fig.  1  Boundary  Configurations 


^  (x,y)  Is  known.  The  density  p  and  the  specific  heat  £  of  the  lubricant  are 
also  assumed  to  be  constants  In  fi 2  •  Furthermore  we  assume  that  the  diffusion  term 
is  much  smaller  than  the  convection  term  along  the  x-axis,  and  that  the  convection 
terra  Is  much  smaller  than  the  diffusion  r 'rm  in  the  y  axis.  From  these  assumptions  the 
original  thermal  energy  equation,  of  an  elliptic  type  in  the  lubricant,  can  be  paraboized 
as  follows.  Let  the  temperature  T; (1=1,2, 3)  be  defined  in  each  corresponding  subdomain 
i) ;  .  Then  the  thermal  energy  equation  in  each  subdcmain  can  be  written  as 


(2) 


—  K|  VZ  Xi  =  0  in  Q;  ,  (i  -  1  and  3) 

-k  i||+fg  =  $<*.»  i”02 


where  y  =  cpu,  and  and  M  is  the  dynamic  viscosity.  As  mentioned 

before  we  assumed  that  the  term  5-y  ]  =  o(l)  and  |CfV  *  o(l) 

in  the  thermal  energy  equation  in  f2a  ,  where  u  and  v  are,  the  velocity  components 
along  the  x-axis  and  y-axis,  respectively. 

Since  the  thermal  energy  equation  in  is  parabolized,  we  impose  the 

boundary  condition  only  on  the  upstream  boundary,  rau  ,  and  proper  juncture 
conditions  on  J(  and  J 1,  i.e.. 


T2=  T  on  Tzv  (3) 

o 

where  the  temperature  of  incoming  lubricant,  T  ,  is  specified.  At  the  juncture 

o 

surfaces  J(  and  J^.  we  require  continuity  of  the  temperature  and  its  normal  heat 
flux,  i.e., 


T  j  -  Tj  and  k |  T  ^  ~  k  j  T ^  on  J 

and  Tz  =  T5  and  kz  T^y  =  k5  T-jy  on  J 


(A) 


In  the  first  numerical  model  problem,  the  following  three  standard  types  of 
boundary  conditions  on  p,  and  are  treated  specifically  for  testing  numerical 
convergence : 

(i)  Dirichlet  type;  T  -  T  on  p-  (i  =  1  and  3)  (5) 

i  o 

where  Te  is  specified. 
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(II)  Neumann  type;  ft.  IT 1  —£..(*  y)  on  F|  ( 1  ■  1  and  3)  (6) 

where  f  (i=l  and  3)  Is  specified, 
i 

(III)  Robin  type;  |31  +  hi  T\  =  h|  T0  on  fi  (i  =  l  and  3)  (7) 

where  the  heat  transfer  film  coefficient  hj  (1=1  and  3)  and  the  ambient  oil  temperature 

distribution  T  (x,y)  along  P(  and  p^  are  specified.  Here  we  assume  the  boundaries, 
o 

P  and  p3  are  in  an  oil  bath  with  an  ambient  oil  temperature  T0(x,y).  It  should 
be  noted  here  that  the  boundary  condition  on  the  downstream  boundary  (outlet)  of 
lubricant,  Pap  >  Is  not  specified  but  is  obtained  as  a  part  of  the  solution.  This 
is  because  the  oiiginal  elliptic  equation  has  been  reduced  to  a  parabolic  equation. 

For  the  purpose  of  the  error  and  convergence  test  in  the  first  model  problem,  we 
take  H]  =  11^  =  H3  =  2 ,  k  \  /  kz  =  3,  k  t  =  k  j  =  3,  h(  =  h  5  =  1,  and  the  juncture  boundaries 
J(  and  Jzare  y  ■  +  1,  respectively.  We  begin  with  the  exact  solution  given  by  a  simple 
polynomial  function  in  each  subdomain  as  follows; 


T,  (x,  y)  - 

_  z 

-2X  +  X 

+  -  3^  +  3 

in  0.| 

T2  = 

-ix*  +  X 

+  Y1  +  Y 

in 

(8) 

T*  Cx,  y )  = 

-lx1  +  x 

+ir+-£r 

in  Cl3 

From  a  given  arbitrarily  specified  function  of  Y  ,  one  can  easily  compute  the  heat 
generation  5.  by  Eqs  (2)  and  (8)  as 

<£  =  YC4-X-1 )  +2  (9) 
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The  three  different  types  of  boundary  conditions  in  this  model  problem 
are  computed  from  the  known  exact  solutions  given  in  Eq(8).  For  example,  with  the  Robin 
type  condition,  hj  =1  (i=l,3),  T0(x,y)  was  computed  from  Eq  (8)  as 

on  p* .  ,  ( i  =  1  and  3) 

hi.  U  x 

and  (10) 

t„  =  +  y  °" 

In  the  second  model  problem,  we  only  treat  the  Robin  boundary  condition 
on  [~ j  and  to  simulate  a  more  realisitic  experimental  condition.  For  this  case, 
the  geometry  of  the  bearing  pad  and  collar  and  the  lubricant  film  thickness  are 
chosen  as  an  analogous  model  in  two  dimensions  corresponding  to  the  three-dimensional 
experimental  condition. 

GALERKIN'S  METHOD  AND  NUMERICAL  PROCEDURE 
Before  we  describe  Galerkin's  finite  element  method  applied  to  the  model 
problem  formulated  In  the  previous  section,  it  is  convenient  to  introduce  a  single 
continuous  temperature  function  T(x,y),  defined  in  the  entire  domain,  Q  =fl|UnzU  SI3  , 
as  follows; 

T(x,y)  =  Tj  (x,y)  in  O;  (1  =  1, 2, 3).  (11) 

To  construct  the  bilinear  functional  in  weak  form,  we  first  introduce  the  test 
function  T*  in  the  test  function  space,  and  next  define  the  inner  product  of  the 
original  partial  differential  equations  In  (2)  and  the  test  function  T*.  By 
integrating  by  parts  the  inner  product  reduces  to 
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fik‘7T'7T* +  |  **  f;  ^  Jn  K3  VT7^cjy<iy 

■^UTtV+J^.  H3TT^  j[$T^ 

+  W,T0T*ds  -v  f  K^T.rJs 
Jr,  'a 


where  the  trial  function  T*  =  0  on  P„,,  and  the  trial  function  T  is  chosen  so  that 
the  essential  condition  T  =  To  on  PiXJ  is  satisfied.  Eq  (12)  is  a  weak  form  for  the 
Robin  condition.  For  the  Dirichlet  type  condition  on  pj  and  p  ,  the  line  integrals 
along  P(  and  P3  in  F.q(12)  are  not  present.  On  the  other  hand  the  trial  function  should 
satisfy  the  Dirichlet  conditions  and  the  test  function  is  chosen  to  be  zero  on  p  and  p  . 
In  the  case  of  the  Neumann  type  condition,  the  boundary  integrals  along  p  and  jp 
appearing  on  the  left-hand  side  of  F,q  (12)  should  not  be  present  and  h  To  in  the 
Integrands  of  the  boundary  integrals  along  p  and  P3  ,  on  the  right-hand  side  of  F.q 
(12),  should  be  replaced  by  f](x,y),  (i  =  l  and  1). 

It  is  of  interest  to  note  that  the  juncture  conditions  on  (J,  and  given 
in  Eq  (4)  are  satisfied  as  natural  conditions  in  Galerkin's  functional  equation 
gi/en  above.  In  the  numerical  computations  an  isoparametric  linear  element  is  used 
as  the  basis  for  both  trial  and  test  functions  throughout  the  present  computations. 

This  choice  of  basis  function  is  equivalent  to  the  centered  finite  difference 
approximation.  In  a  straightforward  manner,  the  bilinear  form  in  Eq  (12)  is  reduced 
to  a  set  of  algebraic  equations.  The  coefficient  matrix  obtained  is  not  symmetric 
but  still  has  a  banded  structure.  The  asymmetry  is  due  to  the  presence  of  the 
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convective  term  In  a  subdomain  .  The  Gaussian  elimination  method  Is  used  to  solve 
the  reduced  matrix  equation. 

NUMERICAL  CONVERGENCE  TEST  IN  THE  FIRST  MODEL 

An  extensive  numerical  test  of  the  convergence  has  been  made  In  the  first 

model  problem.  To  test  numerical  convergence  of  the  present  numerical  scheme, 

we  define  the  error,  E  (1=1,4),  In  four  different  ways  using  the  known  exact 

i 

solution  T  given  in  Eq  (8)  and  the  finite  element  numerical  solution  T(x,y)  as 
follows : 

E,  =  Ht-TIL 

E*=  \  +^K[7tr-T)]l<i>J>y/i 

E,  =  ||vfr-T)||2 

E*=  |lT-tll2 

where  ||  loo  and  |1  jl ^  are  the  well-known  max  norm  and  Lj  norm,  and  defined  as 

Ily-Tll  =  max  IT  ~j\ 

111  ,U°°  ^  (14) 

llT-TlU  = 

and  where  k  =  kj  in  Q.j  (1=1, 2, 3) 

In  the  finite  element  mesh  subdivisions,  Ax/  Ay  =1  is  used  throughout 


the  first  model  problem,  where  Ax,  and  Ay  are  the  lengths  of  the  finite  element 
along  the  x-  and  y-  axes  respectively  (l.e.,  a  square  element  is  used  for  this  model 
problem).  Ten  different  sizes  of  uniform  square  elements  are  tested  in  the  range  of 


I  is  L/A>  <  20  .  The  specific  mesh-subdivisions  tested  are  given  in 

Table  1.  In  the  present  model  problem,  the  computations  are  made  for  three  values 

of  Peclet  number:  Pe  =*  Y  Hz/2  =  0.01,  1,  and  100. 


Table  1  Ten  different  finite  element  subdivisions  used  for  the  error  test. 


Test  case  (i) 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

I 

2 

4 

6 

8 

10 

12 

14 

16 

18 

20 

J 

3 

6 

9 

12 

15 

18 

21 

24 

27 

30 

EL 

6 

24 

54 

96 

150 

216 

294 

384 

486 

600 

N 

12 

35 

70 

117 

176 

247 

330 

425 

532 

651 

MB 

0 

17 

23 

29 

35 

41 

47 

53 

59 

65 

DIM 

108 

595 

1610 

3395 

6160 

10127 

15510 

22525 

31388 

42315 

I  =  Number  of  element  along  the  x-axis 

J  =»  Number  of  element  along  the  y-axis 

EL  =*  Total  number  of  elements 
N  *  Total  number  of  nodes 

MB  =  Bandwidth  (asymmetric) 

DIM  »  Core  memory  space  for  the  coefficient  matrix 


If  we  assume  that  the  error  behaves  like  £  ■  Ax)1"1 

the  limit  Ax— *0,  then  we  may  represent  the  error  as 

Ei  =  A  C'/i)", 

where  A  is  constant  and  I  is  defined  in  Table  1.  By  taking  the  log  of  Eq  (15),  we 
obta I n 

In  Ei  —  In  A  -  ^  In  I 

(16) 

From  our  numerical  results  for  ten  different  finite-element  mesh  subdivisions. 


(  i=l ,2,3,4)  as 

(15) 
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we  have  plotted  the  curve  of  the  points  (  -In  A  ,  In  E |  )  9hown  In  Figs.  2 

through  4.  From  the  two  finest  finite  element  subdivisions,  the  values  of 

index  n  and  the  constant  A  are  obtained  for  three  different  values  of  the  Peclet 

number,  Pe,  and  also  for  three  types  of  boundary  conditions  on  f*i  and  •  The  results 

are  given  in  Table  2  and  3.  In  Table  2  the  values  of  n  for  Ea  and  Ej  are  almost 

-I 

one,  l.e.,  the  corresponding  convergence  error  is  linear,  as  a  function  of  (I)  for 
all  three  types  of  boundary  conditions  and  all  three  Peclet  numbers  tested  except 
for  Pe  =100  in  the  Neumann  boundary  condition.  It  is  surprising  to  see  that  the 
the  convergence  of  the  E  j  error  is  accelerated  as  the  Peclet  number  increases 
-  this  is  contrary  to  our  expectation.  It  is  also  difficult  to  draw  any  conclusion 
on  the  behavior  on  the  convergence  of  the  error  E^.  as  the  Peclet  number  increases. 

In  Table  3,  the  constant  A  for  Et  increases  as  the  Peclet  number  increases  with 
all  three  types  of  boundary  conditions  on  and  .  However,  the  constant  A 

for  E2  and  E3  does  not  vary  much  for  all  three  types  of  boundary  conditions  and 
for  the  three  values  of  Peclet  numbers. 

From  the  results  shown  in  Figure  3  and  Table  2  and  3,  the  present  numerical 
evidence  shows  uniform  convergence  of  the  present  numerical  scheme.  However,  a 
rigorous  mathematical  error  analysis  and  convergence  proof  is  still  open  for  the  class 
of  parabolic-elliptic  coupled  problem  treated  here. 
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Figure  2  (Continued) 


r 


Table  2.  The  Values  of  n  In  various  Error  Norms  for  Three 

Peelet  Numbers,  P  . 

e 


Boundary 

P 

— 

E. 

E 

E 

E, 

Condition 

e 

1 

2 

3 

4 

0.01 

1.0038 

1 .0004 

1.0000 

2.0000 

Robin 

1 

1  .0004 

1.0000 

1.9971 

100 

3.2114 

■HFm 

■tVwwVU 

0.01 

1.2871 

1  .0004 

1.0000 

1.9985 

Neumann 

1 

2.0098 

1.0004 

1.0000 

1.9834 

100 

3.6481 

1.3138 

1.3709 

4.1813 

0.01 

.9175 

1.0004 

1.0000 

2.0000 

Dirichlet 

- 

1 

2.0366 

1.0003 

.  9999 

2.0003 

100 

2.1011 

1.0009 

1.0008 

1.9957 

Table  3.  The  Values  of  Constant  A  in  various  Error  Norms  for 

for  Three  Peelet  Numbers,  P 

’  e 


Boundary 

Condi t ion 

P 

e 

F-1 

E2 

F3 

E4 

Robin 

0.01 

1 

100 

sum 

.  1 721 E+02 
.  1 726E+02 

.  1 60 1 E  +  05 

Neumann 

0.01 

1 

100 

.8260E-01 
.  5876E  +  Q1 
. 1 299E+05 

.4761 E+02 
.4769E+02 
.1281 E  +  03 

PjUiijBWI 

. 1 740E+02 
.  1641 E+02 
. 1257E+06 

Dirichlet 

0.01 

1 

100 

.  1  740E-01 
.6339E+01 
. 1 472E+02 

L - 

.4761 E+02 
.  4768E+02 
.  4823E+02 

.  2993E+02 
.3004E+02 
.  3085E+02 

.  1720E+02 
.  1 747E  +  02 
. 1 8o6E +02 
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RESULTS  FOR  A  SECOND  MODEL  PROBLEM 

For  the  second  model  problem,  the  following  geometrical  and  material  data 
are  used: 

H| =  H3=  0.75  Inch 

H2=  0.001  and  0.0001  inch 

L  =  2.5  inch 

u0=  9.55  inch/sec  and  47.75  inch/sec 
C  =  0.5  Btu/lbm/ eF 
k ,  =  k3  =  26  Btu/hr/ft/°F 
k2  =  0.075  Btu/hr/ f t/°F 
T0  =  100  °F 

h  ,  =  h3  =  30  Btu/hr/f  t2/  °F 
P  =  6.5  X  10  lbm. sec/in 
P  =  0.84  X  10  lbf  sec  /in 
J  =  9336  in.lbf/Btu 

where  U0  *.s  the  maximum  velocity  in  the  lubricant  film.  The  value  of  is  approxi¬ 

mated  by 

5.  <"•''>  ■  “  <  fj  *  * 

and  Y is  computed  by  using  the  mean  velocity,  i.e.,  u  =  ,  since  the 

velocity  is  assumed  linear  between  zero  on  J(  and  up  on 

If  all  the  dimensional  quantities  are  converted  to  consistent  dimensions 
using  (Btu,  sec,  in,°F  ),  then  the  following  values  of  coefficients  are  obtained 
for  use  in  Eqs  (2)  and  (7): 
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ki  =  k  3  =  601  xio 


-  6 


J-Kj 


Ki  =  1-74  *  to 


-6 


J+u 


i*»*  eyrri 


h,  —  hj  =  f7.8?  x  10 


-6 


3+u 


Ssc.  in*  #p 

Table  4  shows  four  specific  test  conditions,  which  are  the  combinations  of  two 
two  velocities  and  two  film  thickness. 


Table  4.  Four  tested  cases  for  the  values  of  Y and  ^  . 


Case  No. 

(inches) 

uQ  (inch/sec) 

y (Btu/sec/in  /°F) 

2£(Btu/in  /sec) 

1 

0.001 

9.550 

0.0810 

0.06427 

2 

0.001 

47.750 

0.0405 

1.6070 

3 

0.0001 

9.550 

0.0810 

6.4270 

4 

0.0001 

47.750 

0.0405 

160.7000 

In  the  present  computations,  two  different  sets  of  mesh  subdivisions  are  used; 
the  first  set  of  data  with  a  coarse  mesh  has  77  nodes  and  60  elements.  The  second 
set  of  data  with  a  fine  mesh  has  315  nodes  and  280  elements.  In  the  fine  mesh, 
we  took  four  uniform  rectangular  elements  In  the  lubricant  (  0,2.)  and  five  uniform 

rectangular  elements  In  both  bearing  pad  and  collar  (  Q ,  and  )  along  the  y-axis, 

and  twenty  elements  along  the  x-axis. 

The  agreement  between  the  computed  temperatures  obtained  by  a  fine  mesh  and  a 
coarse  mesh  was  good.  Therefore  only  the  results  obtained  by  using  the  fine  mesh 
are  shown  in  Fig  5  through  8.  In  this  three  dimensional  computer  plot  of  the  temperature 
distribution,  a  total  of  315  nodes  were  used  with  linear  interpolation.  It  should 
be  noted  here  that  the  film  thicknesses  (i.e.,  0.001  and  0.0001  inches)  were  stretched 
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Figure  5  Temperature  distributions  for  Case  * 
H2~  0.001  inch,  u^  9.55  in/sec 


0.001  inch,  uD=  47.75  in/sec 


Figure  7.  Temperature  distributions  for  Case  3 
H2=  0.0001  inch,  Uq=  9.55  in/sec 


0.0001  inch,  47.75  in/sec 


very  much  along  the  y-axls  for  better  Illustration  in  the  computer  plots. 

The  results  for  all  four  cases  show  that  the  effect  of  the  inlet  temperature 
o 

(100  F  is  used  here)  is  limited  to  a  small  local  region  around  the  inlet.  This 
result  shows  that  there  exists  a  very  thin  thermal  boundary  layer  at  the  inlet  (i.e., 
the  initial  layer  in  the  time-dependent  problem).  This  results  is  not  surprising,  since 
the  heat  is  generated  only  in  the  lubricant. 

It  is  of  interest  to  note  that  for  all  four  cases,  the  temperature  T  at  the 
second  node  on  the  boundary  J|  (or  J*  )  from  the  inlet  (  fiu »  x  =  0  )  Is  lower  than 
that  at  the  adjacent  node  inside  the  solids  (i.e.,  fl(or  Qj  ).  This  means  that 
there  is  a  small  region  of  local  backflow  of  heat  flux.  In  other  words,  in  this 
region,  the  heat  flux  vector  is  pointing  from  the  point  in  the  solid  to  the  lubricant, 
even  though  the  only  heat  generating  source  in  this  problem  is  in  the  lubricant  region. 

CONCLUDING  REMARKS 

From  the  numerical  results  presented  for  the  first  model  lubrication  problem, 
we  can  conclude  that  the  seemingly-unstable  classical  Galerkin  method  is  uniformly 
stable  over  the  range  of  the  Peclet  number  from  0.01  to  100.  This  stability  is 
problably  a  result  of  the  subdomain  of  the  parabolic  equation  being  sandwiched  between 
two  adjacent  subdomains  which  are  elliptic  without  a  convective  term  in  the  heat  trans¬ 
fer  equations.  It  appears  that  the  elliptic  type  equations  for  the  top  and  bottom 
subdomains  play  a  role  in  stablizing  the  numerical  scheme  even  though  we  use  the 
classical  Galerk  n  method  which  is  equivalent  to  the  centered  finite  difference 
approximation.  Ii  the  second  model  problem,  a  local  backflow  of  the  heat  flux  and 
the  presence  of  the  thermal  boundary  layer  are  illustrated.  Future  work  should 
include  the  convective  term  in  the  bearing  collar  where  numerical  stability  may 
not  result  using  the  present  method,  unless  a  proper  weighting  function  is  Introduced 
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in  the  Inner  product. 
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DTNSRDC  ISSUES  THREE  TYPES  OF  REPORTS 

1.  DTNSRDC  REPORTS.  A  FORMAL  SERIES.  CONTAIN  INFORMATION  OF  PERMANENT  TECH 
NICAL  VALUE.  THEY  CARRY  A  CONSECUTIVE  NUMERICAL  IDENTIFICATION  REGARDLESS  OF 
THEIR  CLASSIFICATION  OR  THE  ORIGINATING  DEPARTMENT. 

2.  DEPARTMENTAL  REPORTS.  A  SEMIFORMAL  SERIES.  CONTAIN  INFORMATION  OF  A  PRELIM 
INARY.  TEMPORARY.  OR  PROPRIETARY  NATURE  OR  OF  LIMITED  INTEREST  OR  SIGNIFICANCE 
THEY  CARRY  A  DEPARTMENTAL  ALPHANUMERICAL  IDENTIFICATION. 

3.  TECHNICAL  MEMORANDA,  AN  INFORMAL  SERIES,  CONTAIN  TECHNICAL  DOCUMENTATION 
OF  LIMITED  USE  AND  INTEREST.  THEY  ARE  PRIMARILY  WORKING  PAPERS  INTENDED  FOR  IN 
TERNAL  USE.  THEY  CARRY  AN  IDENTIFYING  NUMBER  WHICH  INDICATES  THEIR  TYPE  AND  THE 
NUMERICAL  CODE  OF  THE  ORIGINATING  DEPARTMENT  ANY  DISTRIBUTION  OUTSIDE  DTNSRDC 
MUST  BE  APPROVED  BY  THE  HEAD  OF  THE  ORIGINATING  DEPARTMENT  ON  A  CASE  BY  CASE 
BASIS. 


